miRNA profiling of chicken follicles during follicular development

MicroRNAs (miRNAs) play a crucial role as transcription regulators in various aspects of follicular development, including steroidogenesis, ovulation, apoptosis, and gene regulation in poultry. However, there is a paucity of studies examining the specific impact of miRNAs on ovarian granulosa cells (GCs) across multiple grades in laying hens. Consequently, this study aims to investigate the roles of miRNAs in chicken GCs. By constructing miRNA expression profiles of GCs at 10 different time points, encompassing 4 pre-hierarchical, 5 preovulatory, and 1 postovulatory follicles stage, we identified highly expressed miRNAs involved in GC differentiation (miR-148a-3p, miR-143-3p), apoptosis (let7 family, miR-363-3p, miR-30c-5p, etc.), and autophagy (miR-128-3p, miR-21-5p). Furthermore, we discovered 48 developmentally dynamic miRNAs (DDMs) that target 295 dynamic differentially expressed genes (DDGs) associated with follicular development and selection (such as oocyte meiosis, progesterone-mediated oocyte maturation, Wnt signaling pathway, TGF-β signaling pathway) as well as follicular regression (including autophagy and cellular senescence). These findings contribute to a more comprehensive understanding of the intricate mechanisms underlying follicle recruitment, selection, and degeneration, aiming to enhance poultry’s reproductive capacity.


Small RNA sequencing and data analysis
Total RNA was extracted from all collected samples using RNAiso Plus reagent (TaKaRa, Japan, #9108), following the manufacturer's instructions.The integrity and quality of the total RNA were assessed using a Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, USA) and an RNA 6000 Nano kit.Three small RNA libraries were constructed for each stage using the QIAseq miRNA Library Kit, following the manufacturer's protocols.The prepared libraries were sequenced on a HiSeq 2500 platform, with a single-end sequencing length of 50 bp, at Novogene Bioinformatics Technology Co., Ltd (Beijing, China).
To obtain high-quality reads, cutadapt-1.2.1 was used to filter low-quality reads and remove reads with adaptor contamination.The remaining high-quality reads were then mapped to the chicken reference genome (GRCg6a) using Bowtie.The reads that did not align to rRNA, tRNA, snRNA, scRNA, and snoRNA were aligned to miRbase to identify known miRNAs in Gallus gallus.The expression levels of miRs were estimated using TPM (transcripts per million).The normalization and identification of differentially expressed miRNAs (DEMs) were calculated using the EdgeR package in RStudio.DEMs with |log FC| > 1 and FDR < 0.05 were considered significant.Furthermore, maSigPro 17 was used to identify miRNAs with significant temporal expression changes across the 10 time points developmentally dynamic miRNAs (DDMs) were defined as DEMs with a goodnessof-fit (R 2 ) of at least 0.5.

Target prediction and functional enrichment
The target prediction and construction of the miRNA-mRNA interaction network construction were detailed in our previous study.In brief, we employed two algorithms, TargetScan (http:// www.targe tscan.org/ vert_ 61/) and miRDB (http:// mirdb.org/) 18 , to predict the potential target genes of the identified miRNAs exhibiting significant expression profile differences over time.Furthermore, we investigated the targeted regulatory relationships between miRNAs and their target differentially expressed genes by retaining miRNAs whose targets overlapped with the DEGs 15 .This analysis was conducted using Cytoscape 19 .For the biological pathway enrichment analysis of DEGs in the four identified clusters, we utilized g:Profiler (https:// biit.cs.ut.ee/ gprofi ler/ gost?tdsou rcetag= s_ pcqq_ aiomsg) with default parameters.The organism option selected was Gallus gallus, and the analysis was run as a multi-query.KEGG enrichment analyses of DDGs for the top expressed miRNAs, and visualization of pathways were performed using OmicShare tools, a free online platform for data analysis (https:// www.omics hare.com/ tools).

Morphological characterization during chicken folliculogenesis
Folliculogenesis refers to the progressive growth and development of ovarian follicles, which includes pre-hierarchical and hierarchical follicles leading up to ovulation or the development of atretic follicles.In our study, we classified the follicular hierarchy based on the color, number, and diameter of ovarian follicles during the peak laying period.Notably, there were significant differences in follicle diameter between adjacent developmental stages (Fig. 1a).During the pre-hierarchical stage, numerous small white follicles (SWF) and large white follicles (LWF) with diameters less than 5 mm were observed, along with 5-6 small yellow follicles (SYF) and large yellow follicles (LYF) ranging from approximately 5-12 mm in diameter.These follicles were categorized as pre-hierarchical follicles.The hierarchical stage consisted of a series of yellow yolky follicles.These preovulatory follicles were characterized by size, with diameters exceeding 20mm.They were denoted as F5 to F1, representing an order of development.During the process of follicle selection, the SYF gradually developed into preovulatory ovarian follicles, progressing from F5 to F1 until ovulation occurred.Additionally, one to two regressing POFs could be observed in the subsequent days. www.nature.com/scientificreports/

miRNA expression profiling patterns during chicken folliculogenesis
We performed a temporal small RNA analysis on chicken GCs collected at 10 stages during folliculogenesis, comprising 4 pre-hierarchical follicles, 5 preovulatory follicles, and one postovulatory follicle.A total of 413.27 million reads (M) of high-quality reads were obtained from thirty small RNA sequencing libraries, with an average depth of 13.78 M per sample.The average GC content across all samples was 48.21%.Furthermore, at least 99.12% of the reads met or exceeded Q20 (Table S1), indicating the high quality of the sequencing data.The high repeatability of each biological replicate was also confirmed by the high Pearson coefficient (Pearson's R > 0.97) within samples per stage (Fig. 1b).Additionally, we observed stronger correlations within groups compared to between groups (Fig. 1c).
Hierarchical clustering analysis revealed consistent patterns between miRNA expression and morphological characteristics across the 10 developmental stages of follicle selection.The progressive phases of ovarian follicles were generally classified into two clusters: pre-hierarchical follicles (enclosed within the larger dotted box) and hierarchical follicles (enclosed within the second largest dotted box).The POFs are distinct from the other stages (enclosed within the smallest dotted box) (Fig. 1d).Interestingly, in our study, the LYFs (also referred to as F6,                                                                                                                                                      www.nature.com/scientificreports/ the smallest preovulatory follicle) appeared to be mixed with F5.We identified 834 known chicken miRNAs expressed in at least one of the samples by aligning against the miRBase database.Among them, approximately 32% were expressed across all 10 stages, with the majority (598) expressed in F3, followed by F2, SWF, and F5 (Fig. 1e).

Expression dynamics of miRNAs throughout chicken folliculogenesis
We observed a total of 176 unique DEMs through pairwise comparisons at each time point.Comparisons between non-consecutive time points resulted in more DEMs than consecutive time points (Fig. 3a, Table S2).Among all stages, POF exhibited the highest number of DEMs (554), followed by 397 DEMs for SWF and 394 DEMs for SYF.Additionally, there were 53 core DEMs that exhibited differential expression across all time points (Fig. 3a,b).The number of DEMs within the pre-hierarchical follicle (40 DEMs), preovulatory follicle (52 DEMs), and post-preovulatory follicle (116 DEMs) periods was relatively low.Notably, the comparison between POF and SWF revealed the smallest number of DEMs, consistent with the HCL tree and Pearson correlation results.Moreover, the number of DEMs increased during the pre-hierarchical follicle period and decreased during the preovulatory follicle period, particularly among the up-regulated DEMs (Fig. 3c,d).

Integrative analysis of DDMs and DDGs during chicken follicle development
To identify the potential targets of the 55 DDMs that changed with time-series during chicken GCs development, we employed a combination of TargetScan and miRDB to predict intersected targets.These predicted targets were then overlapped with 3669 DDGs identified in our previous study using MasigPro 15 .This is just overlapping and a simple correlation, and we did not have any evidence that these miRNAs and genes were negatively or positively correlated.Subsequently, we constructed regulatory networks between miRNA:mRNA pairs using Cytoscape.This analysis revealed significant regulatory relationships between 48 DDMs and 295 target DDGs (Table S4).Notably, gga-miR-363-3p exhibited the largest number of targets (46 DDGs), followed by gga-miR-206-55 (36 DDGs) and gga-miR-106-5p (35 DDGs).Among these 295 targets, 105 DDGs were targeted by 7 of the most highly expressed miRNAs (Fig. 4a).Specifically, gga-miR-148a-3p, which was the most highly expressed miRNA across all 10 stages, targeted 25 DDGs, including PTEN, ING2, RNF38, XPO4, IGF2BP3, AGFG1, and others.

Key signaling pathways during chicken follicle development
To gain insights into the potential biological functions of target genes regulated by differentially expressed miRNAs with dynamic expression patterns over time, we performed functional enrichment analysis.The results revealed the enrichment of pathways associated with various stages of follicle development.For Cluster 1 (decreasing trend), which is characterized by decreasing expression over time, the target genes were enriched in pathways related to follicular development and ovulation, including the Wnt signaling pathway, TGF-β signaling pathway, Oocyte meiosis, and Progesterone-mediated oocyte maturation.Additionally, pathways associated with follicular regression, such as Autophagy and Cellular senescence, were enriched.Cluster 2 (POF dominating trend), characterized by a dominant expression in the postovulatory follicle, showed enrichment in pathways such as Tight junction and ECM-receptor interaction, which are involved in follicular development.Cluster 4 (bell-shaped trend), exhibiting an initial increase and subsequent decrease in expression, showed enrichment in the Ubiquitin mediated proteolysis pathway.Moreover, the target genes of Cluster 3, which showed an increasing trend in expression, were associated with the MAPK signaling pathway (Fig. 4b and Table S5).Furthermore, the target genes of the top 10 highly expressed miRNAs were also enriched in similar pathways.These included Endocytosis, Ubiquitin mediated proteolysis, Notch signaling pathway, and the pathways above related to follicular development, ovulation, and regression (Fig. 4c and Table S6).The network construction of miRNA-mRNA-pathway revealed the enrichment of genes involved in oocyte meiosis, with 31 genes (corresponding to 31 miRNAs), and progesterone-mediated oocyte maturation, with 22 genes (including 28 miRNAs) (Table S7).These pathways are closely related to follicular development and maturation.Among the identified genes, CPEB2/3/4, targeted by 8 miRNAs including gga-let7j-3p, gga-miR-202-5p, gga-miR-363-3p, and gga-miR-92-3p, are known transcription factors involved in oogenesis and spermatogenesis.CDC27, targeted by gga-miR-6700-3p and gga-miR-24-3p, may play a role in regulating the timing of mitosis (Fig. 5).Notably, miR-128-3p, which has been implicated in regulating chicken GC functions through the 14-3-3β/FoxO and PPAR-γ/LPL pathways 47 , was also found to participate in oocyte meiosis (targeting: YWHAB, PPP1CC, and MAPK14) and progesterone-mediated oocyte maturation (targeting: MAPK14 and PDE3B) (Fig. 5 and Table S7).The decreased expression of miR-20b-5p and miR-363-3p may be implicated in the occurrence and progression of human polycystic ovary syndrome (PCOS) 48 .
Pathways related to cellular adhesions play a crucial role in cell motility, proliferation, and differentiation.We identified enrichment of pathways such as Focal adhesion (40 miRNAs targeting 58 mRNAs), Adherens junction (24 miRNAs targeting 24 mRNAs), Tight junction (30 miRNAs targeting 33 mRNAs), and ECM-receptor interaction (24 miRNAs targeting 22 mRNAs), all of which are associated with follicular development (Table S7).

Discussion
Understanding the molecular mechanisms underlying egg-laying performance in hens is crucial for improving poultry breeding strategies.In this study, we performed small RNA-seq analysis on chicken GCs collected from pre-hierarchical, preovulatory, and postovulatory follicles across 10 stages of the laying cycle.We identified a total of 834 known miRNAs, among which 17 highly expressed miRNAs associated with GC development and apoptosis accounted for a significant proportion.Differential expression analysis identified 55 developmentally dynamic miRNAs (DDMs) out of the 176 DEMs.Among these DDMs, 48 miRNAs showed significant regulatory relationships with 295 targeted differentially expressed genes.We gained insights into the functional pathways involving the identified DDMs by constructing co-expression networks of miRNA-mRNA pairs.We predicted potential miRNA-target interactions in pathways related to follicular development, selection, maturation, ovulation, and atresia.Interestingly, many miRNA-mRNA pairs we identified have been previously validated and studied in other experimental settings.
Notably, miR-128-3p, which has been implicated in regulating chicken GC functions through the 14-3-3β/ FoxO and PPAR-γ/LPL pathways 47 , was also found to participate in oocyte meiosis (targeting: YWHAB, PPP1CC, and MAPK14) and progesterone-mediated oocyte maturation (targeting: MAPK14 and PDE3B) (Fig. 5 and Table S5).The decreased expression of miR-20b-5p and miR-363-3p may be implicated in the occurrence and progression of human polycystic ovary syndrome (PCOS) 48 .miR-202-5p in geese has been shown to inhibit lipid deposition and progesterone synthesis by targeting ACSL3 in hierarchical GCs 49 .It also regulates the proliferation and apoptosis of GCs during follicular selection by targeting BTBD10 through the PI3KCB/AKT1 signaling pathway 50 .In chicken, miR-148a-3p upregulation in GCs facilitates steroid hormone biosynthesis by targeting OSBP11 20 .miR-31-5p and miR-202-5p target matrix metalloproteinase and lipid-metabolism-related genes to regulate ovarian development in hens 3 , while miR-20b-5p is associated with the proliferation of chicken www.nature.com/scientificreports/primordial germ cells (PGCs) 51 .miR-26a-5p promotes theca cell proliferation by targeting TNRC6A 29 , and miR-99a inhibits cell proliferation by targeting SMARCA5 in disease-infected chickens 43 .miR-143 is involved in the regulation of cell proliferation and apoptosis 33 , it negatively regulates steroid hormone synthesis and secretion, and its inhibition promotes GC differentiation and follicle development by binding to FSHR 24 .The miR-30 family plays a critical role in cell proliferation and is upregulated in broody chicken GCs.miR-30a-5p can alleviate follicle atrophy by inhibiting autophagy and apoptosis of chicken GCs by targeting Beclin1 20 .miR-101-3p binds to the immune-related gene IRF4 to decrease pro-inflammatory cytokines 41 and potentially modulates ovarian development by targeting BMP5 11 .miR-449b-5p regulates steroidogenesis and estrogen secretion by targeting IGF2BP3 52 , and circadian miR-449c-5p targets ATPZB4 to regulate the calcium signaling pathway 53 .The target genes of miR-7 are enriched in key reproductive pathways that may be linked to egg-laying traits, as they modulate primordial follicle activation, growth, and ovulation 54 .The miR-132 family is implicated in affecting the sexual maturity of laying hens by regulating lipid metabolism through targeting related genes 55 .miR-1662 affects egg laying by regulating metabolism and immunity 56 .miR-128 is associated with lipid regulation during follicular selection in geese 12 .
Our functional enrichment analysis identified several follicular development-related pathways targeted by multiple miRNA-mRNA-pathway networks.For instance, Oocyte meiosis and Progesterone-mediated oocyte maturation have been implicated in the regulation of follicular development and maturation 3,15 .Autophagy and Cellular senescence pathways have been associated with follicular atresia 7,15 .The TGF-beta and Wnt signaling pathways have been found to be involved in follicular selection 11,57,58 .TGF-β is believed to regulate the expression of Hsd3b during ovulation, which is crucial for progesterone (P4) synthesis.P4, derived from steroid hormones, plays a vital role in avian reproduction 57 .The GCs of the most recently ruptured follicle (POF) exhibit a robust www.nature.com/scientificreports/secretory capacity for P4 7 , and POF may serve as a transient supplementary endocrine gland in the chicken ovary, stimulating the development of the pre-hierarchical follicles 59 .Due to the POFs can promote the proliferation of theca externa cells of SWF 59 , we observed an exclusive phenomenon that the number of up-regulated miRNAs are higher than down-regulated miRNAs in POF and SWF comparison.However, in our study, EMC receptor interaction, which is essential for maintaining normal ovarian function and follicle development 60 , GnRH secretion, and Cholesterol metabolism involved in follicular maturation 49,58 , and the Calcium signaling pathway implicated in ovulation 61 were not significantly enriched.Multiple regulatory networks exhibit crosstalk at various stages, involving the targets of gga-miR-202-5p, gga-miR-20b-5p, gga-miR-128-3p, and gga-miR-21-5p.These miRNAs are involved in pathways such as Oocyte meiosis, Progesterone-mediated oocyte maturation, Cellular senescence, Autophagy, TGF-beta signaling pathway, and cellular junction, which play roles in follicular development, maturation, selection, and atresia.miR-202-5p is crucial for follicular development and selection as it inhibits lipid metabolism and steroidogenesis 62 and promotes apoptosis while suppressing the proliferation of geese GCs 50 .Its significant upregulation in large follicles also influences the final maturation of goat-dominant follicles 63,64 .gga-miR-20b-5p is associated with the proliferation of chicken primordial germ cells (PGCs) 51 .Moreover, it is involved in lipid metabolism, and its overexpression can alleviate adipocyte differentiation in the context of human PCOS 48,65 .gga-miR-128-3p regulates follicular selection by promoting GC apoptosis, inhibiting lipid synthesis, and reducing the secretion of progesterone and estrogen 47 .In rat GCs, miR-128-3p improves follicular development by inducing apoptosis 32 .The expression of gga-miR-21-5p increases with the follicle maturation, peaking in POF.miR-21 has been found to have anti-apoptotic effects 66 and contributes to the involution of sheep atretic follicles 67 .It is positively associated with oocyte maturation, and the knockdown of miR-21 in GCs increases apoptosis and is associated with reduced ovulation rate 10,[21][22][23]34 .

Conclusions
We conducted miRNA profiling of GCs from peak laying hens across 10 grades, representing pre-hierarchical, preovulatory, and postovulatory follicles stages.Through this analysis, we identified a substantial number of DDMs and DDGs associated with follicular development, maturation, selection, and regression.These results provide a solid foundation for further investigations into the underlying mechanisms of reproductive performance, shedding light on global functional genes and essential pathways involved in domestic poultry.Additionally, these findings serve as a valuable resource for poultry breeding, aiming to enhance egg-laying capacity.

Figure 1 .
Figure 1.Expression profiles of miRNAs across 10 grades of follicles.(a) Morphological characterizations of ovarian follicles at different stages.(b) The Pearson correlation coefficient shows the repeatability within each sample.(c) Pearson correlation coefficient comparing the grades of follicles (pre-hierarchical, preovulatory, and postovulatory follicle) within and between groups.(d) Hierarchical clustering analysis of miRNA expression patterns.(e) Statistics on the number of miRNAs expressed at each stage.

Figure 3 .
Figure 3. Expression dynamics of miRNAs during chicken folliculogenesis.(a) Comparison of DEMs through pairwise comparison at each time point.(b) Venn diagram showing the interactive dynamics of DEMs at each time point, with the number in the middle representing the core DEMs expressed across all time points and the numbers around the petals indicating the total number of DEMs expressed at each time point.(c) Histogram illustrating the distribution of up-and down-regulated DEMs in POF.(d) Venn diagram displaying the interactive dynamics of DEMs between different grades of follicles (pre-hierarchical, preovulatory, and postovulatory follicles), with the number in the middle representing the core DEMs expressed across the three grades and the numbers around the petals indicating the number of unique DEMs for each grade.(e) Heatmap showing the normalized expression levels (Z-scores of TPM) for each DDM across the four clusters.(f) Temporal expression profiles of the four clusters, with grey lines representing expression levels for individual DDMs, and blue lines representing the mean expression levels for all DDMs within the respective cluster during folliculogenesis.

Figure 4 .
Figure 4. Target genes and pathways of developmentally dynamic miRNAs.(a) Sankey diagram showing the target genes of the top 7 highly expressed miRNAs.(b) Significantly enriched pathways associated with the target genes of all dynamic differentially expressed miRNAs.(c) Pathways enriched by the target genes of the top 10 highly expressed miRNAs.

CPEB2Figure 5 .
Figure 5. Schematic representation of miRNA-mRNA pairs and pathways related to follicular development.

Figure 6 .
Figure 6.Schematic representation of miRNA-mRNA pairs and pathways involved in the follicular selection.

Figure 7 .
Figure 7. Schematic representation of miRNA-mRNA pairs and pathways associated with follicular atresia.

Table 1 .
The top 10 most highly expressed miRNAs, as well as their biological functions.TPM represents the mean expression level of a single miRNA across 30 libraries, the Rate of TPM represents the proportion of a single miRNA's mean expression level compared to all miRNAs.